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ANALYSIS AND MITIGATION OF NUMERICAL DISSIPATION IN 
INVISCID AND VISCID COMPUTATION OF VORTEX-DOMINATED FLOWS 

Osama A. Kandil 

Accomplishments 

In the period of February 1988 to November 1988, The Principal Investigator with the 
assistance of one of his Ph.D. Students has achieved the following accomplishments: 

1. Publication: Kandil, O.A. and Chuang, H. A., “Unsteady Ddta-Wing How 

Computation Using An Implicit Factored Euler Scheme,” AIAA 88-3649-CP, 
AIAA/ASME/SIAM/APS 1st National Fluid Dynamics Congress, Cincinnati, Ohio, July 
25-28, 1988, pp. 248-255. A copy of the paper is attached. 

Abstract 

The conservative unsteady Euler equations for the flow relative motion in the moving frame 
of reference are used to solve for the steady and unsteady flows around sharp-edged delta wings. 
The resulting equations are solved by using an implicit approximately-factored finite-volume 
scheme. Implicit second-order and explicit second- and fourth-order dissipations are added to 
the scheme. The boundary conditions are explicitly satisfied. The grid is generated by locally 
using a modified Joukowski transformation in cross-flow planes at the grid chord stations. The 
computational applications cover a steady flow around a delta wing whose results serve as the 
initial conditions for the unsteady flow around a pitching delta wing about a large angle of attac . 
The steady results are compared with the experimental data and the periodic solution is achieved 
within the third cycle of oscillation. 

2. Journal Publication: Kandil, O. A. and Chuang, H. A., “Computation of Vortex- 
Dominated Flow for a Delta Wing Undergoing Pitching Oscillation,” AIAA Journal, 
Vol. 28, No. 9, September 1990. A copy of the manuscript is attached. 

3. Publication: Wong, T. C., Kandil, O. A. and Liu, C. H., “Navier-Stokes Computations 
of Separated Vortical Flows Past Prolate Spheroid at Incidence,” AIAA Paper No. 
89-0553, Reno, Nevada, January 9-12, 1990. 


Abstract 

The problem of steady incompressible viscous flow past prolate spheroids at incidence is 
formulated using the unsteady incompressible thin-layer Navier-Stokes (NS) equations and the 
unsteady compressible thin-layer Navier-Stokes equations. For the incompressible equations, a 


certain level of unsteady artificial compressibility is added to the continuity equation to secure 
the coupling with the momentum equations during the psuedo-time stepping. The two sets of 
Navier-Stokes equations are solved using a psuedo-time stepping of the implicit flux-difference 
splitting scheme on a curvilinear grid, which is generated by a transfinite grid generator. The 
Baldwin and Lomax algebraic eddy viscosity model is used to model the turbulent flow. The 
computational applications cover a 6:1 prolate spheroid at different angles of attack and different 
Reynolds number. The results are compared with the experimental data. 
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Abstract 

The conservative unsteady Euler equations 
for the flow relative motion in the moving frame 
of reference are used to solve for the steady and 
unsteady flows around sharp-edged delta wings. 

The resulting equations are solved by using an 
Implicit approximately-factored fimte-vo ume 
scheme. Implicit second-order and explicit 
second- and fourth-order dissipations are added 
to the scheme. The boundary conditions are 
explicitly satisfied. The grid is generated by 
locally using a modified Joukowski transformation 
in cross flow planes at the grid chord stations. 
The computational applications cover a * e as 

;ir "nti %s;?ir 9 rihr^^i 

of ° attack Pl t The " 9 teady* rSultHre th 

e experimental data" and the periodic solution 
is achieved within the third cycle of oscilla- 
tion. 

Introduction 

-Unsteady flows around delta wings are 
characterized by the existence of unsteady large- 
and small-scale vortices 

nnssible tretiary vortices), moving shock waves 
Sith time-dependent strengths, time-dependent 
vortex-core formation and breakdown (bubble and 
spiral -vortex breakdown), and interaction of 
Shock waves with the vortical -region and the 
surface-boundary-layer flows. ^! e wit h 9 the 

wi ng 63 ^ ructu ral^Vesponse Tau^ing aeroelastic 

■ ss^opS 11 \r^rrth5>ss ri 5 

C e^tic ti 0 an a llysis tea h d a y s " VeTTsLssed , “.nil 

recoLendations for future code development for 
separated and vortex-dominated flows 
presented. 

The literature on the computational solution 

and experimental data of the unstea y t ranS onic 
dominated flows, particularly i the transon c 
reqime, is unfortunately very limited. ™\ s 
attributed to the complexity of the flow and its 
dependence on numerous parameters, and the 
substantial computational cost involved for the 
flo*T resolution and the time-accurate computa- 
tions . 
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Most of the existing unsteady computational 
schemes are based on the unsteady small distur- 
bance (UTSO) theory 2 * 4 , unsteady full potential 
( UFP) equation 5 * 7 , UTSD.equation with non-isen- 
tropic flow corrections 8 and UFP equation wit 
non-isentropic flow corrections . These schemes 
are restricted to attached flows only. For 
mildly separated flows, integral and fln 'be- 
difference boundary-layer schemes have been 
coupled with potential flow schemes * • 

The unsteady Euler equations adequately 
model shock waves and their motion entropy 
increase across shocks and entropy gradient and 
vorticity production and convection behi d 
shocks, as can be seen from Crocco's theorem and 
the inviscid vorticity transport equation. More- 
over the computational solution of Euler equa- 
tions adequately models separated flow from s arp 
oH( .., 12-14 For smooth-surface separation, 

round-edge' separation, shock-induced separation, 

viscous diffusion and dissipation, vortex break- 
down flow transition and turbulence; viscou 
terms must be added to Euler equations to recover 
the full Navier- Stokes equations or an approxi- 

for. of the*. editions. « tt.ogh o 
process of adding the viscous term 
unsteady Euler solvers is simple and straight 
forward the computational cost for g 
Reynolds-number unsteady flows is substantia 
might be prohibitive due to the need for tine 
grids to adequately resolve the viscous effects. 

Euler/Navier-Stokes zonal approaches 1 ’ 
have been demonstrated to maintain the accuracy 
of the Navier-Stokes solutions and in the mean 
time alleviate to a good extent the computa- 
Mnnil rost of the Navier-Stokes equations. 

rsrssa. "ss. a* £a. r ir “ 

further for unsteady flows. 

Recently successful time accurate solutions 
of the unsteady Euler and Navier-Stokes ^ua- 
have been presented for airfoils > • 

The only existing unsteady Euler solutions for 
vorteT-domi nated flows are those of the rol ing- 
, r j nation of a sharp-edge delta wing in a 
locally conical supersonic flow around a mean 
anqle of attack and a zero angle of attack, which 
were presented by the authors in Refs. 13 and 14. 
The authors derived the unsteady Euler equations 
for the flow relative motion in a moving frame of 
reference, and the equations have been solve i y 
uslna an explicit, multi-stage time stepping, 
f n 9 te-volumn P scheme. Periodic solutions were 
achieved in the third cycle of rolling oscilla- 
H Details of the surface pressure, cross- 

flow* velocity and cross-flow Mach contours were 
presented C show1 ng the primary vortex and wave 
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shocks formation and interaction. 

In the present paper, the three-dimensional 
unsteady Euler equations in a moving frame of 
reference are solved by using an implicit 
approximatel y-f actored finite-volume scheme. The 
unsteady airfoil computations of this scheme were 
compared earlier with the experimental data of 
ref. 21, which were in good agreement. The 
present three-dimensional steady results are 
compared with the experimental data of ref. 22 
for two levels of explicit dissipation. With the 
steady results serving as initial conditions for 
the unsteady flow, the flow around the same delta 
wing undergoing pitching oscillation about the 
quarter-chord axis is solved in this paper. 

Formulation 
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Starting with the conservative form of the 
unsteady Euler equations for the flow absolute 
motion in the space-fixed frame of reference, and 
using the following relations for the substantial 
and local 
vector "A" 
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Cartesian coordinates x', y' , z' of the moving- 
frame of reference, the resulting Euler equations 
are 
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relative motion; E , F and G the inviscid 
r r r 

fluxes of the relative motion, $ a source term 
due to the motion of the reference frame, p the 
density, p the pressure, e and h the total energy 
and total enthalpy per unit mass, V and a the 
flow absolute velocity and absolute acceleration, 
V and a p the flow relative velocity and rela- 
tive acceleration, V and a^ the transformat ion 

velocity and transformation acceleration, V and 
- o 

the translation velocity and translation 

acceleration of the moving frame, Hi and Hi the 
angular velocity and angular acceleration of the 
moving frame, r the position vector of a fluid 
particle with respect to the moving frame, and 
y is the ratio of specific heats. The pitch, yaw 
and roll angles are referred to by using the 
Eulerian angles a, 6 and e; respecti vley . The 
" 1 " refers to the derivatives, time or coordi- 
nates with respect to the moving frame of 
reference. 

Introducing the curvilinear coordinates £ 1 , 
n’ and c' in the moving frame of reference, which 
are given by 


Z* s Z' (x‘ ,y' ,z' ) , n' 8 n'(x' ,y‘ ) » 
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Eqs. (2)-(7) are transformed to 
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In Eqs. (20) and (21), A r , B r> C r and H r are the 

Jacobian matrices aE p / aq f , aF p / aq p , aG p / aq p and 

aS/aq r , respectively; 5 . , 6 ^ . and 6^. the three- 

point central difference operators; D , , D m ^, 

and D . the implicit dissipation operators; and 
me r 

D , , D ■ and D . the explicit dissipation 
operators. The expressions of the dissipation 
operators are given in ref. 24. The implicit 
damping coefficient is and the explicit 

damping coefficients are and The damping 

coefficients are having the same values in the 
, n' and c 1 di rections . 


pw r W r + i' g , P, o w r h r l t 


s = 


The solution of Eq. (20) is obtained through 
(18) three successive sweeps in the n 1 * c‘ and^ 5 ' I 

respectively. Once Aq 1 ^ is obtained, q f is 
m q i found from 


In Eqs. ( 15) -( 19) . J’ 1 is the Jacobian of 
transformation and U r , V r and W r are the contra- 
variant velocity components. 


Conputatlonal Scheme 

Equation (14) is integrated over 5 ', n' and 
£* the divergence theorem is applied to the 
resulting equation. By using the implicit 
approximate factorization scheme, we obtain the 
following difference equation for a typical cell 
(i , J * k) 
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The surface boundary conditions is explicitly 
enforced through the normal momentum equation 24 
while the farfield boundary conditions for 
subsonic flows are enforced using the inflow- 
outflow conditions which are based on the Rieman 

i nvariants 2 *. 


Computational Results 

A sharp-edged delta wing of aspect ratio, 
AR , of one, at a mean angle of attack, o^, of 
20.5° and in a free stream Mach number, of 

0.3 is considered for the computational 
application of the implicit three-dimensional 
vectorized program. The body conformed grid 
consists of 80x38x48 cells in the s' , V and 
s' directions, respectively; and its size is one 
root-chord ahead of the wing vertex, two root- 
chords behind the trailing edge and one root- 
chord radius In the cross flow planes. The outer 
boundary consists of a hemi-spherical surface 
with its center at the wing vertex and a cylin- 
drical surface with its axis coinciding with the 
wing axis. The grid is generated in cross-flow 
planes using a modified Joukowski transformation 
which is locally applied at the grid chord sta- 
tions with exponential clustering at the wing 
surface. 

Steady Flow: 

The Implicit program is used, to solve for 
the steady flow at 20.5° angle of attack with two 
levels of numerical dissipation; a low dissipa- 
tion (L0) with e 2 = 0-05. e 4 - 00.0025 and 
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e m = 0.25 and a high dissipation (HD) with = 

0.25, = 0.0025 and = 0.25. Figure 1 

shows the solutions at two chord stations of 0.52 
and 0.81 on the wing. The figures froin left to 
right show the surface pressure coefficient 
(Figures 1.1), the static pressure coefficient 
contours (Figures 1.2) and the cross-flow 
velocity directions (Figures 1.3). Here, we 
compare the surface pressure of the Implicit- 
Scheme with low dissipation with that of the 
experimental data of Hummel^. Other comparisons 
are given by the authors in ref. 24. 

At X'/C = 0.52, the computed surface 

pressure with low dissipation is in excellent 
agreement with the experimental data. The 
computed static pressure contours with low dissi- 
pation show higher pressue levels than those 
computed contours with high dissipation, parti- 
cularly in the vortical core region. With low 
dissipation, the highest pressure contour is 1.6, 
while with high dissipation the highest pressure 
contour (at the same locatin) is 1.1. Comparison 
of the computed surface pressure with high dissi- 
pation with that of the experimental data (given 
in ref. 24) shows that the computed peak suction 
pressure is underpredicted by 13% although the 
remainder of the surface pressure is in good 
agreement with the experimental data. 

At X'/C = 0.81, the computed surface pres- 
sure with low dissipation is higher than that of 
the experimental data, particularly under the 
primary vortex core. The computed peak suction 
pressure is about 25% higher than that of the 
experimental data. Comparison of the computed 
static pressure contours with low dissipation 
shows higher pressure levels than those computed 
contours with high dissipation, particularly in 
the vortical core region. The cross flow veloci- 
ties of both dissipation levels (Figures 1.3) 
show almost identical shapes and directions. 

Figure 2 shows the exper imental static 
pressure contours of Hummel^ in planes 
perpend icul ar to the wind direction (Figures 

2.1) , the computed static pressure contours in 
planes perpendicular to the wing surface (Figures 

2.2) and the cross-flow velocity directions 
(Figures 2.3). The computational results along 
with the experimental data are shown for two 
cross flow planes in the wake; X'/C - 1.02 and 
X'/C = 1.25. At X'/C = 1.02, the computed outer 
contours are in excellent agreement with the 
experimental contours. For the most inner static 
pressure contours, the experimental data show 
higher level than those of the computed results. 
On the other hand, the implicit results with low 
dissipation show higher level than those of the 
higher dissipation. Similar results are seen at 
X'/C = 1.25. At this location, it is seen that 
the trai 1 i ng-edge vortex core is captured using 
the low-dissipation implicit scheme. 

The discrepancies between the experimental 
data and the computed results with low dissipa- 
tion level are attributed to the grid coarseness 
in the vortical core and to the viscous effects 
in the vortex core as well as on the wing upper 
surface. The discrepancies between the computed 
resuKs with low and high dissipation are 
obviously due to the low and high values of the 


explicit second-order damping coefficient. 

On the VPS-32 computer of NASA Langley 
Research Center, a typical steady flow case takes 
1050 psuedo time stepping to reach a residual 
error of 10*^. 

Unsteady Flow (Pitching Oscillation about the 
Quarter-Chord Axis): 

The steady results with low dissipation 
level are used as the initial conditions for 
calculating the unsteady flow around the same 
wing which is undergoing a pitching oscillation 
about the quarter chord axis. The angle of 
attack a(t) is given by 


a(t) = a m + a 0 Sin 2/y kt 

where a is the amplitude, and k is the reduced 

0 ★ 

k C * 

frequency (k = k = dimensional frequency 

and c = wing chord“length) . In this application 
~ 20.5°, a Q - 2°, = 0.3 and k = 3 which 

corresponds to a period of 2.95 per cycle. Each 
cycle of oscillation takes about 1,475 time steps 
and the solution covers 5,000 time steps which 
correspond to 3.39 cycles of oscillation. Figure 
3 shows a vz t motion at the top, which is 
followed by the surface pressure variation, 
static pressure contours and cross-flow velocity 
in each row of figures. The numbers 1-15 on the 
a-t curve and on the other figures show the 
instants at which the computational results are 
shown. Here, we show the computations at X'/C = 
0.52 and the computed surface pressures are shown 
every 200 time steps starting from the 2,200 time 
step, which corresponds to point 1 on the a-t 
curve. The static pressure contours and the 
cross-flow velocity directions are given at the 
3,000; 4,000 and 5,000 time steps which corre- 
spond to points 5, 10 and 15, respectively; on 
the a-t curve. Comparison of the surface 
pressure at points 7 and 14, corresponding to 
2.31 and 3.25, cycles respectively, shows that 
periodic oscillation has already been reached. 
Experimental data for unsteady vortex dominated 
wing flows are urgently needed for bench-mark 
conpari sons . 

Concluding Remarks 

The three-dimensional unsteady Euler 
equations in a moving frame of reference are 
solved by using an implicit approximately- 
factored finite-volume scheme. The computational 
applications cover steady low-subsonic flow 
around a sharp-edged delta wing at a large angle 
of attack, and unsteady low-subsonic flow around 
the same wing undergoing a pitching oscillation 
about the same large angle of attack. The steady 
flow problem has been computed with two levels of 
numerical dissipation, and the results have been 
compared with each other and with the experi- 
mental data. The low-dissipation results give 
better agreement with the experimental data than 
those of high-dissipation results. However, fine 
grid embedding are needed in the vortical 
regions. Moreover, the viscous terms must be 
added in the vortical regions (free-shear layers) 
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and in the surface boundary-1 ayer flow. Outside 
of these regions, Euler equations are adequate 
for computing the flow field. This calls for 
urgent extension of the Euler/Navier-Stokes zonal 
scheme 16 to this problem. The results of the 
unsteady flow application show consistency, and 
they show periodic solution in the third cycle of 
oscillation. Although the code has been verified 
previously with unsteady airfoil computation 1 ^, 
unsteady experimental data for vortex dominated 
flows are urgently needed for bench-mark compari- 
son. Work is underway to use the flux-difference 
splitting scheme with the unsteady Euler/Navier- 
Stokes zonal scheme. 
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Figurel . Comparison of the surface pressure*stat1c pressure and cross-flow 
velocity at two chord stations on the wing for a sharp-edge delta 
wing, 80x38x48 cell , M = 0.3, a ■ 20.50®, AR * 1; 1. Surface 

pressure, 2. Static pressure contours, 3. Cross-flow velocity. 
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FIGURE 3. PITCHING MOTION, VARIATION OF SURFACE PRESSURE, STATIC-PRESSURE 
CONTOURS AND CROSS-FLOW VELOCITIES FOR A SHARP-EDGED DELTA WING 
UNDERGOING PITCHING OSCILLATION; 80X38X48 CELL, = 0.3, 

„ = 20.50. a = 20° . k = 3. AR = 1. e„ = 0.05, e, = 0.0025, 






To At?*** ~ 

AIAA ‘3eOV»— f’. 


COMPUTATION OF VORTEX-DOMINATED FLOW FOR A DELTA 
WING UNDERGOING PITCHING OSCILLATION 1 


Osama A. Kandil* and H. Andrew Chuang** 

Old Dominion University, Norfolk, Virginia 

Abstract 

The conservative unsteady Euler equations for the flow relative to a 
moving frame of reference are used to solve for the three-dimensional steady 
and unsteady flows around a sharp-edged delta wing. The resulting equations 
are solved by using an implicit, approximately-factored, finite-volume 
scheme. Implicit second-order and explicit second- and fourth-order 
dissipations are added to the scheme. The boundary conditions are explicitly 
satisfied. The grid is generated by locally using a modified Joukowski 
transformation in cross-flow planes at the grid chord stations. The 
computational applications cover a steady flow around a delta wing whose 
results serve as the initial conditions for the unsteady flow around a 
pitching delta wing at a large mean angle of attack. The steady results are 
compared with the experimental data and the unsteady results are compared with 
results of a flux-difference splitting scheme. 


f Thi s paper has been presented as AIAA-88-3649-CP, Cincinnati, Ohio, 1988. 

*Professor, Department of Mechanical Engineering and Mechanics, Associate 
Fellow AIAA. 
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Introduction 


Unsteady flows around delta wings are characterized by the existence of 
unsteady large- and small-scale vortices (primary, secondary and possible 
tertiary vortices), moving shock waves with time-dependent strengths, time- 
dependent vortex-core formation and breakdown (bubble and spiral types of 
vortex breakdown), and interaction of shock waves with the vortical -region and 
the surface-boundary-layer flows. These highly unsteady aerodynamic loads may 
interact with the wing structural response causing aeroelastic instabilities 
which may restrict the aircraft flight envelope. In Ref. 1, the status of 
computational unsteady aerodynamics for aeroelastic analysis has been 
discussed, and recommendations for future code development for separated and 
vortex-dominated flows are presented. 

The literature on the computational solution and experimental data of the 
unsteady vortex-dominated flows, particularly in the transonic regime, is 
unfortunately very limited. This is attributed to the complexity of the flow 
and its dependence on numerous parameters, and the substantial computational 
cost involved for the flow resolution and the time-accurate computations. 

Most of the existing unsteady computational schemes for airfoil and wing 

2 - 

flow applications are based on the unsteady small disturbance (UTSD) theory 
unsteady full potential (UFP) equation®"*^, UTSO theory with vorticity and 
non-i sentropi c flow corrections® and UFP equation with non-isentropic flow 

corrections^. These schemes are restricted to attached flows only, and hence 
they cannot be used to capture vortex-dominated flows. For mildly separated 
flows, integral and finite-difference boundary-layer schemes have been coupled 
with potential flow schemes 10,11 . However, such schemes cannot be used for 
massive separation. 


On the other hand, the unsteady Euler equations adequately model shock 
waves and their motion, entropy increase across shocks and entropy gradient 
and vorticity production and convection behind shocks, as can be seen from 
Crocco's theorem and the inviscid vorticity transport equation . Moreover, tne 
computational solution of Euler equations adequately models separated flows 
from sharp edges 1 ^ -14 . For smooth-surface separation, round-edge separation, 
shock-induced separation, viscous diffusion and dissipation, vortex breakdown, 
flow transition and turbulence; viscous terms must be added to Euler equations 
to recover the full Navier-Stokes equations or an approximate form of these 
equations. Although the process of adding the viscous terms to the unsteady 
Euler solvers is simple and straightforward, the computational cost for high- 
Reynolds-number unsteady flows is substantial and might be prohibitive due to 
the need for using fine grids to adequately resolve the viscous effects. 

Euler/Navier-Stokes zonal approaches 15, 16 have been demonstrated to 
maintain the accuracy of the Navier-Stokes solutions and in the meantime, 
alleviate to a good extent the computational cost of the Navier-Stokes 
equations. These approaches along with fine grid embedding in the vortical 
regions should be developed further for unsteady flows. 

Recently, successful time accurate solutions of the unsteady Euler and 
Navier-Stokes equations have been presented for airfoi 1 s 14, 1 ^~^. The only 
existing unsteady Euler solutions for vortex-dominated flows are those for the 
rolling-oscillation of a sharp-edge delta wing in a locally conical supersonic 
flow around a mean angle of attack and a zero angle of attack, which were 
presented by the authors in Refs. 13 and 14. The authors derived the unsteady 
Euler equations for the flow relative to a moving frame of reference, and the 
equations have been solved by using an explicit, multi-stage Runge-Kutta time 
stepping, finite-volume scheme. By casting the equations in a moving frame of 
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reference, we eliminated the computation of the grid motion, Moreover, since 
the equations are expressed in terms of the flow vector field relative to the 
moving frame, the equations still preserve the conservation form. Periodic 
solutions were achieved in the third cycle of rolling oscillation. Details of 
the surface pressure, cross-flow velocity and cross-flow Mach contours were 
presented showing the primary- vortex and shock-wave formation, interaction 
and disappearance. 

In the present paper, the three-dimensional unsteady Euler equations in a 
moving frame of reference are solved by using an implicit approximately- 
factored, finite-volume scheme. The resulting code has been validated through 
the solution of unsteady flows around airfoils undergoing pitching oscilla- 
tion 17 , and comparison of the computed results with the experimental data of 
ref. 21. The present three-dimensional steady results are compared witn the 
experimental data of ref. 22 for two levels of explicit dissipation. With tne 
steady results serving as initial conditions for the unsteady flow, the flow 
around the same delta wing undergoing pitching oscillation about the quarter- 
chord axis is solved in this paper. The results of the pitching delta wing 
are compared with those obtained using the three-dimensional flux-difference 
splitting scheme of ref. 20. 

Formulation 

Starting with the conservative form of the unsteady Euler equations for 
the flow relative to a space-fixed frame of reference, and using the following 
relations for the substantial and local derivatives of a scalar "a" and a 
vector "A" 
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we obtain the conservative unsteady Euler equations for the flow relative to 


moving-frame of reference. In terms of the Cartesian coordinates x 
of the moving-frame of reference, the resulting Euler equations are 
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In Eqs. (1)— (12), q r is the flow vector field of relative motion; E^, and 
G the inviscid fluxes of the relative motion, S a source term due to the 
motion of the reference frame, p the density, p the pressure, e and h the 
total energy and total enthalpy per unit mass, V and a the absolute 

velocity and absolute acceleration of the flow, and a p the flow relative 

velocity and relative acceleration, and a t the transformation velocity 

and transformation acceleration. V q and l Q the translation velocity and 

translation acceleration of the moving frame, w and u the angular velocity 
and angular acceleration of the moving frame, r the position vector of a fluid 
particle with respect to the moving frame, and y is the ratio of specific 


heats. 


The pitch, yaw and roll angles are referred to by using the Eulerian 

angles a, 6 and 9; respecti vley . The refers to the derivatives, time or 

coordinates with respect to the moving frame of reference. Equation (2) shows 
that the present formulation is in a conservation form. However, the 
conservation form is not a strong one due to the existence of the source term 
on the right hand side. But on the other hand, the Jacobian correspondi ng to 
this term is not spatially differenced (see Eqs. 21 and 22) and hence the 
truncation error of this term is mainly due to temporal differencing. 

Introducing the curvilinear coordinates t' , h‘ and C‘ in the moving 
frame of reference, which are given by 

l' = S'(x' ,y 1 ,z' ), V = V(x' ,y' ,z‘ ), 

C = C (x 1 ,y ' , z 1 ) (13) 

Eqs. (2)-( 7) are transformed to 
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In Eqs. ( 15 ) - ( 19 ) , J" 1 is the Jacobian of transformation and U r , V r and are 
the contravariant velocity components. 


Computational Scheme 

Equation (14) is integrated over £' , V and C 1 and the divergence 
theorem is applied to the resulting equation. By using the implicit 
approximate factorization scheme, we obtain the following difference equation 
for a typical cell (i, j, k) 
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In Eqs. (20) and (21), A r , B r , C p and H f are the Jacobian matrices f£ r /oq r> 

OF /dq , 5G / 5q and SS/dq , respectively; 5-., 6 , and & r , the three-point 

central difference operators; 0 m ^. , D , and 0 m{ .. the implicit dissipation 

operators; and D , D . and D the explicit dissipation operators. The 
K et, en ec 

expressions of the dissipation operators in the K‘ direction are given by 
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where U r , V r and W r are the contravariant components of the relative velocity. 
The implicit damping coefficient is e m and the explicit damping coefficients 
are and e 4 . The damping coefficients are having the same values in the 
V , V and C directions. 

The solution of Eq. (20) is obtained through three successive sweeps in 
the V , C' and ; respectively. Once Aq^ is obtained, q" + ^ is found from 


-n+1 -n . -n 

q, = q + &L 


( 28 ) 


Id 


Boundary Conditions 

The surface boundary condition is enforced explicitly through the normal 
momentum equation 

ap * ~ 

— = p V • (V • 7 n) - p n • a (29) 

3n r r t 

where n is the unit normal of the wing surface. Also, the farfield boundary 
conditions are enforced explicitly. In the present application, subsonic flow 
in the farfield is considered and hence, the inflow-outflow boundary condi- 
tions are based on the Riemann invariants, and R^, for one-dimensional 
flow normal to the boundary which are given by 
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co op y- I oo co 

2 . , .1/2 

R. = V. • n + — - ( y p./p. ) (31) 

11 y-1 li 

* 

where n is the unit normal of the outer boundary of the computational region 
and the subscripts ® and i refer to the farfield conditions and the values 
extrapolated from the interior cells at the boundary, respectively. Thus, the 
inflow boundary conditions are given by 
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(V • t) = (V • t) 

(32. d) 


b 


where b refers to the boundary and t is a unit vector tangential to the 
boundary. For the out-flow boundary conditions, the subscript ® in Eqs. 
(32. c) and (32. d) is replaced by the subscript i. Equations (32. a) and (32. b) 
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give a complete definition of the flow at subsonic boundaries. Symmetric flow 
conditions are used at the plane of symmetry for the present application. 

Computational Results 

A sharp-edged delta wing of aspect ratio, AR, of one, at a mean angle of 
attack, a m , of 20.5° and in a free stream Mach number, M o , of 0.3 is 
considered for the computational application of the implicit three-dimensional 
vectorized program. The body conformed grid consists of 80x38x48 cells in 
the V , V and C' directions, respectively; and its size is one root-chord 
ahead of the wing vertex, two root-chords behind the trailing edge and one 
root-chord radius in the cross flow planes. The size of the computational 
domain has been determined through numerical experiments, where the maximum 
absolute value of the difference of the pressure coefficient did not exceed 
0.5S from that of a computational domain which has double the present 
computational size. The small computational domain is used because of the 
limited available computational resources. The outer boundary consists of a 
hemi-spherical surface with its center at the wing vertex and a cylindrical 
surface with its axis coinciding with the wing axis. The grid is generated in 
cross-flow planes using a modified Joukowski transformation which is locally 
applied at the grid chord stations with exponential clustering at the wing 

surface . 

Steady Flow: 

The implicit program is used to solve for the steady flow at 20.5 angle 
of attack with two levels of numerical dissipation; a low dissipation (LD) 
with e 2 = 0.05, e 4 = 00.0025 and e m = 0.25 and a high dissipation (HD) 
with e 2 = 0.25, e 4 = 0.0025 and e m = 0.25. As a general principle, the 
dissipation level in a computational scheme must be as low as possible. But 


it must be enough to obtain a stable solution. Although we are dealing with a 
low subsonic low, where shocks do not exist, second-order dissipation terms 
are required due to the large gradients in the vortical region. The purpose 
of this numerical experiment is to determine a low dissipation level, which 
gives a stable solution also. Figures 1 and 2 show the solutions at two chord 
stations of 0.52 and 0.81 on the wing, respectively. The figures from left to 
right show the surface pressure coefficient (Figures 1.1 and 2.1), the static 
pressure coefficient contours (Figures 1.2 and 2.2) and the cross-flow 
velocity directions (Figures 1.3 and 2.3). Here, we compare the surface 
pressure of the Implicit-Scheme with low dissipation with that of the 
experimental data of Hummel 22 . Other comparisons are given by the authors in 
ref. 24. 

In Figure 1, the computed surface pressure with low dissipation is in 
excellent agreement with the experimental data. The computed static pressure 
contours with low dissipation sn:w higher pressure levels than those computed 
contours with high dissipation, particularly in the vo-tical core region. 
With low dissipation, the highest pressure contour is 1.6, while with high 
dissipation the highest pressure contour (at the same location) is 1.1. The 
cross-flow velocities of both dissipation levels show almost identical shapes 
and directions. Comparison of the computed surface pressure with high 
dissipation with that of the experimental data (given in ref. 24) shows that 
the computed peak suction pressure is underpredicted by 18% although the 
remainder of the surface pressure is in good agreement with the experimental 
data. 

In Figure 2, the computed surface pressure with low dissipation is higher 
than that of the experimental data, particularly under the primary vortex 
core. It should be noted here that the experimental data of Hummel shown in 
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Fig. 2.1 is the average of his data at X'/C - 0.7 and 0.9 in ref. 22. The 
computed peak suction pressure is about 25‘i higher than that of the 
experimental data. Comparison of the computed static pressure contours with 
low dissipation shows higher pressure levels than those computed contours with 
high dissipation, particularly in the vortical core region. With low dissipa- 
tion, the highest pressure contour is 1.7, while with high dissipation the 
highest pressure contour is 1.2. The cross flow velocities of both 
dissipation levels show almost identical shapes and directions. 

Figures 3 and 4 shows the experimental static pressure-coefficient 
contours of Hummel 22 in planes perpendicular to the wind direction (Figures 
3.1 and 4.1), the computed static pressure-coefficient contours in planes 
perpendicular to the wing surface (Figures 3.2 and 4.2) and the cross-flow 
velocity directions (Figures A. 3 and 4.$). The computational results along 
with the experimental data are shown for two cross-flow planes in the wake; 
X'/C = 1.02 and X'/C = 1.25. In Figure 3, the computed static pressure 
contours of the outer vortex-core region are in excellent agreement with the 
experimental contours. For the most inner static pressure contours (levels 
higher than 0.3), the experimental data show higher level than those of the 
computed results. On the other hand, the results of static pressure contours 
of the implicit scheme with low dissipation show higher pressure levels than 
those with the high dissipation. Similar results are seen at X’/C = 1.25 

(Figure 4). At this location, it is seen that the outer contour of the 
trai 1 ing-edge vortex core is captured using the low-dissipation implicit 
scheme (contour level of 0.4). The cross-flow velocities at X /C 1.25 also 

show that the trail ing-edge vortex core has been captured. 

The discrepancies between the experimental data and the computed results 
with low dissipation level are attributed to the grid coarseness in the 


vortical core and more important to the viscous effects in the vortex core as 
well as on the wing upper surface. The discrepancies between the computed 
results with low and high dissipation are obviously due to the low and high 
value of the explicit second-order damping coefficient. 

On the VPS-32 computer of NASA Langley Research Center, the CPU time is 
40 m- sec. per grid point per time step and a typical steady flow case takes 
1050 psuedo time stepping to reach a residual error of 10 . 

Unsteady Flow (Pitching Oscillation about the Quarter-Chord Axis): 

The steady results with low dissipation level, best level for a stable 
solution, are used as the initial conditions for calculating the unsteady flow 
around the same wing which is undergoing a pitching oscillation about the 
quarter-chord axis. The angle of attack a(t) is given by 

a(t) = a + a sin M kt (33) 

mo ® 

★ 

k C * „ 

where a Q is the amplitude, and k is the reduced frequency (k = -^g- , k = 

CD 

dimensional frequency and c = wing chord length). In this application = 
20.5°, = 2°, = 0.3 and k = 3 which corresponds to a period of 2.95 per 

cycle. Each cycle of oscillation takes about 1,475 time steps and the solu- 
tion covers 5,000 time steps which correspond to 3.39 cycles of oscillation. 

Figure 5 shows a vz t motion at the top, which is followed by the 
surface pressure variation, static pressure contours and cross-flow velocity 
in each row of figures. The numbers 1-15 on the a-t curve and on the other 
figures indicate the instants at which the computational results are shown. 
Here, we show the computations at the chord station X'/C = 0.52 and tne 
computed surface pressures are shown every 200 time steps starting from the 
2,200 time step, which corresponds to point 1 on the art curve. The static 


pressure contours and the cross-flow velocity directions are given at the 
3,000; 4,000 and 5,000 time steps which correspond to points 5, 10 and 15, 
respectively; on the a-t curve. Comparison of the surface pressure at 
points 7 and 14, corresponding to 2.31 and 3.25, cycles respectively, shows 
that periodic oscillation has already been reached. 

Considering the time-history of the surface pressures, which are 
indicated by the numbers 1-5, 6-10 and 11-15, the variation of the surface 
pressure and the motion of the primary vortex are predicted as the wing 
pitching motion ( a vz t curve) is progressing. For example, the curves 1-5 
show that as the angle of attack increases from instant 1 to instant 2, tne 
peak suction pressure increases and the primary vortex moves outboard in the 
spanwise direction. This is physically expected since the primary vortex 
strength increases. As the angle of attack increases from instant 2 to instant 
3, the peak suction pressures reaches a maximum value and decreases indicating 
that the peak suction pressure is leading the wing motion. As the angle of 

attack decreases from instant 3 to instant 5, the peak suction pressure 

decreases and the primary vortex moves inboard in the spanwise direction. 

On the surface pressure curves, we also show comparisons at selected 

instants with the computations of Rumsey^. Rumsey's computations are 
produced using a three-dimensional, implicit, flux-difference splitting 
program, which is known as "CFL3D" . The results are computed by using our 
grid which is used to compute our results. The results are in good agreement 
and hence the comparison adds confidence to the computational results. 
Experimental data for the unsteady vortex-dominated wing flows are urgently 
needed for bench-mark comparisons. However, it is emphasized here that our 
present scheme and code have been validated earlier^ for unsteady flows 

around pitching airfoils using available experimental data. A 
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Concluding Remarks 

The three-dimensional unsteady Euler equations in a moving frame of 
reference are solved by using an implicit approximately-factored finite-volume 
scheme. The computational applications cover steady low-subsonic flow around 
a sharp-edged delta wing at a large angle of attack, and unsteady low-subsonic 
flow around the same wing undergoing a pitching oscillation about the same 
large angle of attack. The unsteady application presents pioneering results 
for the first time that we know of. The steady flow problem has been computed 
with two levels of numerical dissipation, and the results have been compared 
with each other and with the experimental data. The low-dissipation results 
give better agreement with the experimental data than those of high- 
dissipation results. However, fine grid embedding are needed in the vortical 
regions. Moreover, the viscous terms must be added for accurate computations 
of the vortex-core regions (free-shear layers) and for the surface boundary- 
layer flow computations as well. Outside of these regions, Euler equations 
are adequate for computing the flow field. This calls for urgent extension of 
the Euler/Navier-Stokes zonal scheme" 0 to this problem. The results of the 
unsteady flow application show consistency, and they show periodic solution in 
the third cycle of oscillation. The present unsteady results have been 
compared with those of the flux-difference splitting scheme, and they are in 
good agreement. Although the code has also been verified previously with 
unsteady airfoil computation^, unsteady experimental data for vortex 
dominated flows are urgently needed for bench-mark comparison. 
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Abstract 

The problem of steady incompressible viscous 
flow past prolate spheroids at incidence is 
formulated using the unsteady incompressible 
thin-layer Navier-Stokes (NS) equations and the 
unsteady compressible thin-layer Navier-Stokes 
equations. For the incompressible equations, a 
certain level of unsteady artificial compressi- 
bility is added to the continuity equation to 
secure the coupling with the momentum equations 
during the psuedo-time stepping. The two sets of 
Navier-Stokes equations are solved using a 
psuedo-time stepping of the implicit flux- 
difference splitting scheme on a curvilinear 
grid, which is generated by a transfinite grid 
generator. The Baldwin and Lomax algebraic eddy 
viscosity model is used to model the turbulent 
flow. The computational applications cover a 6:1 
prolate spheroid at different angles of attack 
and different Reynolds number. The results are 
compared with the experimental data. 

Introduction 

The flow about slender bodies over a wide 
range of angles of attack is characterized by the 
presence of large scale vortices on the leeward 
side of the body as a result of the surface flow 
separations. As the angle of attack changes from 
low to very high values, the three-dimensional 
boundary-layer flow changes from an attached and 
vortex-free flow to a separated and vortex- 
dominated flow. At moderate to high angles of 
attack, the flow separates on the leeward side of 
body forming two symmetric vortices, which in 
turn cause secondary and possible tertiary 
vortices. At higher angles of attack, asymmetric 
steady or unsteady vortex flow develops and 
unfavorable changes occur in the force and moment 
characteristics. For additional details, the 
reader is referred to a recent survey paper by 
Newsome and Kandil 1 , which covers the physical 
aspects and the numerical simulation of 

vortical flows past slender bodies and highly 
swept wings over a wide range of angle of attack 
and Mach number. 
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The present problem has been considered 

computationally 2 ' 13 and experimentally 14 ' 1 by 
several investigators due to its significance to 
aircraft and missiles applications. In particu- 
lar, the problem of incompressible flow past 
prolate spheroids over a wide range of angles of 
attack has been considered for computational 

simulation by several investigators » ’ “ * 
due to the existence of detailed experimental 
data 14 ' 17 . Three-dimensional integral boundary- 
layer equations have been used by Stock and Tai 
to solve for the flow around slender porlate 
spheroids at incidence. Ragab 4 has used the 
three-dimensional incompressible boundary-i ayer 
equations on a non-orthogonal coordinate system 
to solve for laminar, transitional and turbulent 
flows past a 6:1 ellipsoid of revolution at 

various incidences. In a recent paper, Ragab 
coupled the boundary-layer equations with an 
inhomogeneous form of Euler equations to effect a 
mechanism of vorticity generation. The method 
was applied to the turbulent flow past a b.i 
prolate spheroid at incidence. Boundary layer- 
equations have also been used by Yanta and 
Wardlaw 3 , Jettmar and Kordulla 3 and Patel and 
8reak 7 . 


Newsome and Adams^ jsed the MacCormack s 
explicit scheme to solve for the vortical flow 
over elliptical body missile using the unsteady 
compressible thin-layer Navier-Stokes equations. 
No skin-friction coefficient calculation or 
comparison was shown. Hartwich and Hall solved 
for the vortical flow over a tangent-ogive 
cylinder using an implicit flux-difference 
splitting scheme for the unsteady incompressible 
thin-layer Navier-Stokes equations. The computer 
program* used is called "V0R3DI" - the same 
program we used in the present paper for the 
incompressible Navier-Stokes equations. Again, 
no skin friction calculation or comparison was 
shown in their paper. 

Pan and Pulliam 10 used the implicit factori- 
zation scheme to solve the compressible thin- 
layer Reynolds-averaged Navier-Stokes equations 
for the incompressible flow over a 6:l-prolate 
spheroid at 10° angle of attack 0.029 Mach 
number and 1.6xl0 6 Reynolds number based on the 
ellipsoid length. The results were compared with 

the experimental data of Meier, et. al ’ ®jd 
the sting-support of the body in the experiments 
was not modeled in the computational simulation. 
Moreover, geometrical singular lines were present 
along the body axis due to the grid used. e 
computations covered a laminar flow and a turb 
lent flow triggered at a quarter of the body 
length from the body nose where Baldwin and Lomax 
algebraic eddy viscosity model was used. 
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thermal conductivity (in the case of compressible 
Nav i er-Stokes equations) k with 

\i = ^ + u- 3 VI + 
e t 


* FT72TT ,i2) c 1 * *2 “ ]t 1111 

In Eqs . (8)-(10), the pressure is related to the 
total energy per unit mass, e, the kinetic energy 
per unit mass and the density by the ideal gas 
equation 


P => (r-l)p [e - j (u 2 + v 2 + w 2 )] (12) 

In Eqs. ( 2 ) - ( 1 1 ) , the contravariant components of 
the velocity; U, V, U, the metric functions ♦, 
and * 2 and the Jacobian of transformation, J, are 
gi venSy 


U = 5 x u + S y v + 5 z w ( 13 ) 


k e = k > k t = 


(1 


I y 

t r , 


( 20 ) 


where P is the effective viscosity, k 0 the 
effective 6 thermal conductivity, p p the turbuient 
viscosity, P r the laminar Prandtl number, P rt the 
turbulent Prandtl number and C p the specific heat 
under constant pressure. The turbulent viscosity 
p is obtained by using the two-layer algebraic 
eddy viscosity model which was first developed by 
Cebeci^ for the boundary -1 ayer equations and 
modified later by Baldwin and Lomax 20 for the 
Navier-Stokes equations. The turbulent viscosity 
in this model is given by 


v = n x u + n y v + n z w 

(14) 

W = C x u + C y v + : z w 

(15) 

-2 . r 2 r 2 
*1 = ‘’x + S + C Z 

(16) 

♦z 1 («, f S v c + c z w c )/3 

(17) 


J' 1 = x £ (y n i. - y c z„) + x n (y c z r) 
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V>0 


r < r. 


r > r 
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( 21 ) 


where r is the normal distance from the body 
surface (for the prolate spheroid it is in the 
direction of a ray from the axis) and r z the 


smallest value of r at 
turbulent viscosity (n 


Vi 


which the inner-layer 
is equal to the outer- 

1 ayer turbulent viscosity (^ t ) Q * For the inner 
layer, the turbulent viscosity is calculated by 
usina the Van Driest algebraic formula 


+ Z h " y h 


(18) 




P 


.22) 


In Eqs. ( 2 ) - ( 5 ) , a certain level of aritificial 
compressibility is added to the continuity equa- 
tion through the time dependent term (p/P) 
with constant B. The purpose of this term is to 
re-couple the continuity equation with the momen- 
turn equation producing a parabolic set of partial 
differential equations in time during the psuedo- 
time marching. Once steady flow solution is 
reached, the artificial compressibility term 
disappars. The incompressible equations are 
nondimensionalized using the reference parameters 
of t, V., L/V^, and for the length, 

velocity, time, density and molecular viscosity, 
respectively. For the compressible equations, 
the reference parameters are L, a a , L/a^, 
and p for the length, velocity, time, density 
and moTecul ar viscosity, respectively. For the 
two sets of equations, the Reynolds number is 
defined as R e 3 

Turbulence Modeling: 

For turbulent modeling, the Navier-Stokes 
equations are transformed to the Reynolds- 
averaged equations by replacing the coefficient 
of molecular viscosity, n, and the coefficient of 


where the mixing length * is given by 

* ■ k r [1 - exp (- r*/ A*] (23) 


where k is the von Karman constant, A is a 
damping constant and r + is given by 


r * = r (p w V 1/2/p w 


(24) 


The subscript w refers to the body surface. In 
Eq. (22), |wj is modulus of vorticity which is 

given by 

2 >2 

|u| - |^xV| = [(u - v x ) + (v 2 - w y ) 


* (w x - u 2 r] 


y 

2il/2 


(25) 
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>i- 1/2 ]/^ (35) 


and 6- and v are defined similarly. The 
invisci J d fluxes E n , F" and G n at the interfaces 
i * 1/2, j 1 1/2 and k ? 1/2 separate the nega- 
tive and positive waves and hence they can be 
written as the solution to a Riemann problem in 
the form 




“i +1/2 ' i 


(36) 


Similar expressions are used for 1/2 * F j ±1 /2 
and G? l/z . Following Roe's scheme, we rewrite 
(4 E- ±1/ / [similarly for (A F jtl/2 ) n and 
(<’l /2 ) n ] as 

K ±1/2 ) n * ± ( A iti/2 )n s iti/2 qn (37) 


Q? + Q" , 

,.+ *n « + n r v i i + ll 
(A i + l/2 } - A 1 2 ' 


(38) 


Using a similarity transformation, the matrices 
A, B and C are split. For matrix A, the split- 
ting is given by 

A = R A L = R [A + - A*] L * A + - A' (39) 

where ^(‘t U| )/2 are diagonal matrices of 
eigenvalues, R and L are the right and left 
eigenvector matrices of A. Using Eqs . lab) 
(39) into Eq. (34), we obtain 

[I/J i,j,k At) + (A Vl/2 6 i-l/2 


[I/(j i,j,k 4t) + (A i-l/2 6 i -1/2 
‘ A i +1/2 A i+l/2^ n AQ 
s {R i+l/Z A i +1/2 ^ * °' ? ' 5 i + l 
* L i+l/2 } " 6 i+l/2 ^ 

- {R i _i/2 *1-1/2 [I + 0-5 ( *i 

- «;.,)] 6 i.uz ? 

where 

3 diag ($*, * 4 ^ 

and 
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The viscous flux G v is centrally differenced. 
Thus, for the thin-layer Navier-Stokes equations, 
the difference equations becomes 


- A l+l/2 6 i+l/2 ) + (B j-l/2 6 j-l/2 
* B i+l/2 6 i+l/2 ) * (C k-l/2 6 k-l/2 

- w w» B Aqn - - [{A i-i/2 6 i-i/2 

- A i +1/2 *i+l/2 ) + (B j-l/2 6 j-l/2 


CI/(J i,j,k At) + A i-l/2 6 i -1/2 ' A i”+l/2 6 i+l/2 
+ 8 j -1/2 6 j-l/2 * B j '*■1/2 6 j-l/2 
+ ( C+ + Z ^ k -1/2 6 k-l/2 

- (O' * Z) k+1/2 fi k+l/2^ n Aq " - * RES «"> (43) 


‘ B j+1/2 6 j+l/2^ + ^ C k -1/2 6 k -1/2 
‘ C k+l/2 6 k+l/2 )3 Q 


Since the scheme of Eq. (40) is highly 
dissipative, a high resolution scheme is used to 
enhance the scheme accuracy. For one-dimensional 
equation of Eq. (40), the high resolution scheme 
gives 


nere Z is the Jacobian matrix corresponding to 
he viscous flux G y , and RES (Q n ) is the discrete 
epresentation of the spatial derivatives in Eq. 
1) evaluated at n with the high resolution 

cheme applied to the inviscid fluxes. Equation 
431 is solved by using a symmetric, planar 
auss-Seidel relaxation in the 5-direction and 
pproximate factorization in the n and direc- 
ions. This factors Eqs (43) in the form 
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windward angle range ® = 0 - 50°. This is 
attributed to the fully turbulent flow assumption 
of the comoutations while in the experimental 
case the flow is laminar to transional. Beyond 
this angle, the CFL3D results are in excellent 
agreement with the experimental data. The 
results of the Y0R3DI are unacceptable. Finally, 
we show a comparison of the surface pressure and 
skin friction coefficient in the axial vertical 
plane. The C p results of the two codes are close 
while the Cf results are substantially different. 

The CFL3D code took 3000 time steps to 
reduce the residual error by four orders of 
magnitude while the V0R3DI code took 400 time 
steps to reduce the residual error by two orders 
of magnitude. 

2. Fixed Transitional Flow (X/L ■ 0.2), 

Re - 7.7x10®, ■ = 10°: 

Here, the turbulent model is turned on at 
the axial location X/L = 0.2. Figures 10-14 
fully cover this case. In Fig. 10, we show the 
cross-flow velocity vectors in the cross-flow 
planes which are computed by the CFL30. It is 
seen that the flow separation is very small in 
the forward planes, while in the rearward planes 
it grows larger. Figure 11 shows a comparison of 
the computed surface pressure by CFL30 and V0R3DI 
in the same three cross-flow planes. The results, 
are in good agreement. In Fig. 12, the corre- 
sponding computed skin-friction coefficient by 
CFL30I and V0R30I along with the experimental 
data are presented. Again CFL3D results are in 
good agreement with the experimental data while 
those of the V0R3DI show substantial discrepan- 
cies in the windward side. 

Figure 13 shows a comparison of the computed 
surface pressure and skin friction in the axial 
plane. While the computed surface pressure of 
CFL30 and V0R30I are in good agreement, the 
computed skin friction coefficient show substan- 
tial differences. 

Figure 14 shows a comparison of the conver- 
gence history of CFL30 and V0R3DI. It is seen 
that the CFL30 takes 6000 time steps to reduce 
the residual error by four orders of magnitude 
(which is twice the number of time steps needed 
for the turbulent flow case). The V0R30I still 
takes 400 time steps to reduce the residual error 
by two orders of magnitude. 

3. Laminar Flow Case, Re - 1.6xl0 6 , « *10°: 

The results of this case is covered in Figs. 
15-17. Here, we only consider the results of 
CFL3D which are shown after 10,000 time steps. 
In Fig. 15, we show a comparison of the computed 
skin-friction coefficient with those of the 
experimental data at six cross-flow planes. It 
is seen that the results are in good agreement 
with the experimental data in the 9 range of 0° - 
150 in the forward part of the body and in the 
e range of 0 - 115° in the rear part of the 
body. In the remainder range of of 9, the 
computed results are substantially different from 
the experimental data. This is attributed to the 
flow transition from laminar to turbulent in the 
region of separated flow, which has not been 


accounted for in the computational simulation. 

It should be pointed out here that the skin- 
friction results of this c*se obtained here in 
this paper, after 10,000 iterations, are a little 
different from those obtained by the same code in 
Ref. 11, where 6,000 iterators have been used. 

It is obvious from Fig, 17 that the residual 
error is still coverg^ng vie 10,000 iteration 
step after it showed constant residual error in 
the range of 5,000-6,000 iterations. Figure 16 
shows good agreement of the surface pressure in 
the axial vertical plane with that of 
experimental data. The computed skin friction in 
the leeward side shows substantial difference 
with the experimental data. 

Concluding Remarks 

The problem of steady incompressible viscous 
flow past a 6:l-prolate spheroid at incidence has 
been solved using two sets of the Navier-Stokes 
equations; an incompressible set and a compress- 
ible set. The computational scheme used for the 
two sets is the implicit flux-difference split- 
ting scheme. The incompressible set is solved 
using a computer code known as "V0R3DI" and the 
compressible set is solved using a computer code 
known as H CFL3D" . The compressible code has been 
applied to three test cases; a fully turbulent 
flow case, a fixed transition flow case and a 
laminar flow case, while the incompressible code 
has been applied to the first two cases. The 
predicted results of the compressible code are in 
good agreement with the experimental data while 
those of the incompressible code are unacceptable 
with the exception of the surface pressure. For 
the high-angle-of-attack turbulent flow case, the 
incompressible code predicts small primary and 
secondary vortices which are closer to the body 
with larger azimuthal angles than those of the 
experimental data . 

The compressible code, however, takes large 
number of time steps to converge. The number of 
time steps are inversely proportional to the 
Reynolds number. On the other hand, the incom- 
pressible code takes small number of time steps 
to converge. Additional work is needed to modify 
the incompressible flow code for better predic- 
tion of the viscous effects; e.g., skin-friction 
coefficient, details of the primary vortex core 
and the secondary vortex. Perhaps one should 
switch the upwind flux-difference splitting 
scheme to central differencing as the residua 
drops down to order of 10' 2 , since the 
incompressible equations then are dominantly 
elliptic. Obviously, one also need better 
transitional and turbulent models for this 
problem than the simple minded approach of 
switching on the turbulent model at a prescribed 
position. It should be also recalled that some 
of the constants of the two-layer, algebraic 
model of Balwin and Lomax have been obtained from 
matching its results with transonic flow results 
-- a point which needs further investigation. 
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Figure 7. Surface Pressure in Cross-Flow Planes, a=30.0, Re=7,200,000 
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Figure 11. Surface Pressure in Cross-Flow Planes, a=10.0, Re-7,700,000, 
(0.2 Fixed Transition) 
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Figure 12. Skin-Friction Coefficient in Cross-Flow Planes, a=10.0, Re=7,700,000 
(0.2 Fixed Transition) 



Figure 13. Surface-Pressure and Skin-Friction Coefficients in the Axial Vertical 
Plane, a=10.0, Re=7,200,000 (0.2 Fixed Transition) 
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